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ABSTRACT 

The illumination of an accretion disk around a black hole or neutron star 
by the central compact object or the disk itself often determines its spectrum, 
stability, and dynamics. The transport of radiation within the disk is in general 
a multi-dimensional, non-axisymmetric problem, which is challenging to solve. 
Here, I present a method of decomposing the radiative transfer equation that 
describes absorption, emission, and Compton scattering in an obliquely illumi- 
nated disk into a set of four one-dimensional transfer equations. I show that 
the exact calculation of the ionization balance and radiation heating of the ac- 
cretion disk requires the solution of only one of the one-dimensional equations, 
which can be solved using existing numerical methods. I present a variant of 
the Feautrier method for solving the full set of equations, which accounts for 
the fact that the scattering kernels in the individual transfer equations are not 
forward-backward symmetric. I then apply this method in calculating the albedo 
of a cold, geometrically thin accretion disk. 

Subject headings: accretion - radiative transfer - methods: numerical 



1. INTRODUCTION 

The structure and dynamics of an accretion disk around a neutron star or a black hole 
is often strongly influenced by its irradiation from the compact object or the disk itself. 
In galactic sources, for example, the optical emission (see, e.g., de Jong, van Paradijs, & 
Augusteijn 1996), the stability (see, e.g., van Paradijs 1996; King, Kolb, & Burderi 1996), 
and the warping of the accretion disk (Pringle 1996; Wijers & Pringle 1999) are probably 
determined by X-ray irradiation from the central source. In extragalactic sources, the ob- 
served broad iron emission lines (see Nandra 1997 for a review), the weak or absent Lyman 
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edges (Sincell & Krolik 1997; Antonucci 1999), as well as the disk warping (Pringle 1996) 
may also be caused by X-ray irradiation of the accretion disk. 

Calculating the structure even of a non-illuminated accretion disk is not trivial, be- 
cause of the complicated physical processes involved, such as the non-LTE character of the 
transport of radiation (see, e.g., Hubeny & Hubeny 1998), and assumptions regarding the 
vertical profile of the viscous heating (see, e.g., Laor & Netzer 1989). The problem becomes 
even more difficult when illumination is taken into account, both because additional physical 
process, such as photoionization (see, e.g., Raymond 1993, Ko & Kallman 1994), become 
dominant, and because of the multidimensional character of the problem. 

The transport of radiation through an accretion disk illuminated by a distant source is 
in general two-dimensional and not axisymmetric. The spectrum of the reflected radiation, 
when only Compton scattering is taken into account, has often been solved in terms of 
the distribution of photon escape probabilities from the reflecting medium (e.g., Lightman 
& Rybicki 1980) or using Green's functions (e.g., Poutanen, Nagendra, & Svenson 1996). 
When absorption and emission processes are taken into account, the problem is often solved 
assuming plane-parallel symmetry (e.g., Sincell & Krolik 1997), or employing variants of 
the A-iteration (e.g., Basko, Sunyaev, & Titarchuk 1974; Matt, Fabian, & Ross 1993) and 
Monte-Carlo methods (e.g., George & Fabian 1991). 

Solving directly the radiative transfer equation in two spatial dimensions is challeng- 
ing, especially when the energies of the photons change by Compton scattering and the 
problem becomes non-local both in space and in photon energy (see, e.g., Mihalas, Auer, & 
Mihalas 1978; Auer & Paletou 1994; Dykema, Klein, & Castor 1996; Dullemond & Turolla 
2000; Busche & Hillier 2000 for methods of solving problems in two spatial dimensions). In 
general, a problem that is two-dimensional in coordinate space requires the solution of the 
radiative transfer equation in five dimensions, two in coordinate space and three in photon 
momentum space. In the case of illumination of a geometrically thin accretion disk, how- 
ever, the problem can be simplified significantly. Because the photon mean free paths in the 
radial and azimuthal directions are much smaller than any characteristic length scale in the 
accretion flow, a geometrically thin accretion disk can be decomposed into a finite number of 
plane-parallel, obliquely illuminated slabs. Solving the radiative transfer equation for each 
slab requires only four dimensions, one in coordinate space and three in photon momentum 
space. When the angular dependence of the interaction cross section between photons and 
matter can be expanded into a finite number of Legendre polynomials, the problem can be 
simplified even further (Chandrasekhar 1960). In this case, the radiative transfer equation 
is equivalent to a finite number of equations over one dimension in coordinate space and 
two dimensions in photon momentum space, i.e., to a finite number of radiative transfer 
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problems in one spatial dimension. These equations can be solved in general using existing 
numerical methods. 

In this paper I derive the system of equations (reduced to one spatial dimension) that 
describe, in an obliquely illuminated slab, isotropic absorption and emission, as well as 
Compton scattering to first order in e/m c and T c /m c , where e is the photon energy and T c 
and m e are the electron temperature and rest mass. I then explore a variant of the Feautrier 
method for solving the resulting transfer equations, in which the scattering kernels do not 
have forward-backward symmetry (see Milkey, Shine, & Mihalas 1975). I illustrate the use 
of the derived equations and method of solution by solving simple problems related to the 
albedo of a cold disk. More detailed calculations of the coupled radiation and gas properties 
using the approach described in this paper will be reported elsewhere. 



2. RADIATIVE TRANSFER 

I solve the radiative transfer problem in a plane-parallel slab illuminated by some arbi- 
trary, external source of photons. I describe all physical quantities in this system using an 
orthogonal, Cartesian reference frame, with its z-axis parallel to the finite dimension of the 
slab. I also set c = &b = h — 1, where c is the speed of light, /cb is Boltzman's constant, and 
h is Planck's constant. 

I assume that the electrons in the slab have density n c (z) and temperature T e (z) <C m e , 
where m e is the electron rest mass. I also assume that the illuminating radiation is a parallel 
pencil beam of net flux nF. Neglecting induced processes, the radiative transfer problem 
is linear and therefore the effect on any arbitrary illumination pattern can be computed 
by summing the solutions obtained for each plane-parallel pencil beam. The direction of 
illumination is described by the vector f = sin d x sin (pix + sin 6\ cos 4>\y + cos 6{Z, where 6\ and 
(pi are the directional angles. 

I describe the radiation field in terms of the monochromatic specific intensity I(z, /i, 0, e), 
where z is the distance from the edge of the slab, /i = cos 9 and are the directional angles 
of the propagation vector I, and e is the photon energy. Because I study non-polarized 
radiation, I have suppressed the dependence of the specific intensity on polarization mode. 
I assume that absorption and emission in the slab are isotropic and denote the absorption 
coefficient by x{ z i e ) an d the source function by S(z, e). 

Keeping only terms to first order in e/m e and T e /m e , the radiative transfer equation 
that describes absorption, emission, and Compton scattering can be written as (Pomraning 
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l-VI(z,l,e) = x(z,e)[S(z,e)-I(z,l,e)]-n e a T [l-2—)I(z,l,e) 

Tfl e 



+n e a T 



J n=0 



2n + 1 
4tt 



P n (i-i')S n I(zJ',e) , (1) 



where <tt is the angle-integrated cross section for Thomson scattering, dfl' is the solid-angle 
element around P n (j ■ I') is the Legendre polynomial of order n, and 
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Defining the energy-dependent optical depth as 



dr(e) 



■ 2e , 

1 n e a T + X 

m, 



dz , 



the relative importance of absorption and scattering as 
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and the redistribution function in the scattering Kernel as 

3 

P (i-i') = J2^nPn(i-i') 



n=0 



with zu n = (2n + l)p s S n , the radiative transfer equation (1) becomes 



(2) 

(3) 

(4) 
(5) 

(6) 



d 1 f - - 

Equation (7) is a second-order, intergrodifferential equation in a four-dimensional phase 
space and is equivalent to a system of four equations in a three-dimensional phase space (see 
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Chandrasekhar 1960, §48.1). This is true because the expansion in Legendre polynomials of 
the redistribution function (6) terminates after the first four terms. 

Following Chandrasekhar (1960), I will write the radiative transfer equation (7) in terms 
of the specific intensity I<i(t, p, 0, e) of the diffuse radiation field, i.e., of the photons that 
have interacted with the gas at least once. Expanding the specific intensity of the diffuse 
radiation field as 



J d (r, p, 0, e) = 7 cT( r > A*» e ) cos t m Oi - <t>)\ » 



(8) 



m=0 



where 0i defines the plane of illumination, the transfer equation for the diffuse field becomes 
equivalent to the system of equations (cf. Chandasekhar 1960, §48.1) 

^/ d m (r, p, e) = J d m (r, p, e) - 5™p a S(r, e) - ± £ w™P?(jx) £ P™(pVT(r, p\ e)dp> 



-{2 - 5™) J2 ^n(-l) m+n Pn^)P^) 



Fe~ T/ ^ , (9) 



where 5™ is Kronecker's delta, 
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and P™(p) are the associated Legendre's functions of the first kind defined by 
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where C\{p) = (1 — ^ 2 ) 1 ^ 2 I\{^). The second-order equation is 
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Finally, the third-order equation is 

= m - 1^(1 - /i 2 ) 3/2 p s f\i - » /2 ? /2 ww 

-§^(1- P 2 f /2 {1- iAT' 2 Ps Fe-^ . (15) 

Equations (12)— (15) are four, second-order, partial differential equations. Because they 
describe the evolution of the diffuse radiation field, they can be solved with the following 
boundary conditions, 

J d m (r = 0, n < 0, e) = J d m (r max , n > 0, e) = , (16) 

where r max is the total vertical optical depth of the slab, and 

J d m (r, /i, e = 0) = J d m (r, ^ e - oo) = . (17) 

In practice, the vertical optical depth of an accretion disk can be very large and, therefore, 
the second of boundary conditions (16) can be exchanged with another condition that is 
easier to handle numerically. 

The first two moments of the specific intensity of the diffuse radiation field, which are 
necessary for calculating the energy and momentum exchange between photons and matter, 
can be calculated as 



Mr, e) = \ J' I° d (r, e, 



p)dp (18) 
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and 

H d (r, e) = H^t, e) cos^x + H^t, e) sin(0i)y + H 2 (t, e)z , (19) 

where 

Hi(r,e) = \ Jjl-^I^e^W (20) 

and 

^(t, c) = t Z) / ^(r, e ' ^ • ( 21 ) 

m=0 

3. NUMERICAL METHOD 

In §2, I showed that the specific radiative transfer problem in two spatial dimensions 
described by equation (1) has been reduced to four problems in one spatial dimension each, 
which are easier to solve. However, solving the latter problems still requires special care, 
for a number of reasons. First, the interaction of the illuminating radiation with the disk 
material takes place in the outermost layers of the slab, which are optically thin. As a 
result, the method of solution of the transfer equation must be accurate in the limit of 
low optical depth. Second, in a typical disk-illumination problem, the interaction of the 
illuminating radiation with the disk material is dominated by true-absorption at low photon 
energies but is scattering-dominated at high photon-energies. For this reason, simple A- 
iteration procedures are not adequate and either accelerated iterative procedures (e.g., the 
accelerated A-iteration or the method of variable Eddington factors) or other non-iterative 
procedures (e.g., the Feautrier method) must be employed. However, even the latter methods 
are not directly applicable to the problem studied here because the redistribution integrals 
in the right-hand sides of the four transfer equations (12)— (15) are not forward-backward 
symmetric. In this section I describe a variant of the Feautrier method that has the desired 
properties for solving the four one-dimensional radiative transfer equations (12)-(15). I 
follow in general the procedure outlined by Milkey et al. (1975), pointing out the differences 
that arise from the particular properties of the problem studied here. 

I choose as independent variables the Thomson scattering optical depth 

r es (z) = / n e (z)a T dz , (22) 
Jo 

which is independent of photon energy, as well as the photon energy e, and the direction 
of propagation //. I then write the four one-dimensional radiative transfer equations in the 
general form 

Psfj.^- = I m - S m ' s - 5 m ' a m = 0,..,3, (23) 
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where S m ' s and S m ' a are the symmetric and antisymmetric parts of the source functions and 
scattering integrals in equations (12)-(15), which depend implicitly on the radiation field. 
Defining the Feautrier variables 



w m (r es ,e,/i) = -[/ m (r es ,e,/i) + / m (r cs ,e, -fi)] 



(24) 



and 



(r es ,e,/i) -J m (r cs ,e, -pi)} , (25) 
the transfer equation takes the form of the system of equations (Milkey et al. 1975) 

am 



PsP 



Ps/i 



dr cs 



= v 



u m — <S n 



(26) 
(27) 



Note, that because of the lack of forward-backward symmetry in the scattering kernels, 
equations (26)-(27) cannot be combined into a single second-order equation, as in the usual 
Feautrier method (see, e.g., Mihalas 1978). 

I then discretize all quantities over grid points in the variable r es , N e grid points 
in photon energy, and grid points in the direction of propagation ji. For simplicity, I 
use, e.g., u™- to denote the first Feautrier variable of order m of the diffuse radiation field, 
evaluated on the i— th grid point in the quantity r cs , on the k— th grid point in photon energy, 
and on the /— th grid point in the direction of propagation /i, such that j — k + (l — l)iV M . In 
order to recover the diffusion of photons in energy space because of Compton scattering by 
thermal electrons, I use a second-order differencing scheme in photon energy. For example, I 
denote the first and second derivatives of the first Feautrier variable with respect to energy 
by 

du m 1 

" =J2 X ^<) + k (28) 
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no = "2 
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(34) 
(35) 

In differencing with respect to the variable r es , I note that the quantity is density-like 
and I, therefore, use a center-differencing scheme for all interior grid points, i.e., 



2 



= u i+l,j u i,j 

T i+\ T i 

i+1/2 es 'es 



(36) 



On the other hand, the quantity is flux-like and I therefore use for all interior grid points 



dv" 1 



dr . 



, U i+l/2,j U i-l/2,j 
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Using the above differencing schemes, the difference equations in all interior grid points 
become 

1 
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where I have written explicitly the dependence of the source functions on the radiation field 

as 



«+l/2,i 



U i+l/2,jk V i+l/2,k + U i+l/2, : 



(40) 
(41) 



Note here that the above source functions have a different structure than those of Milkey et 
al. (1975), because of the presence of the antisymmetric term G" +1 / 2 j that does not depend 
on the diffuse radiation field. 

At the illuminated surface of the slab, which I denote by % = 1, the boundary condi- 
tion (16) translates into u 7 ^ = and applying equation (38) on the first half of the first 
grid cell, denoted by (1,3/2), I obtain 



V 3/2,j ~ U lj + 



ljk U lk ^lj 



(42) 
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V N d -l/2,j 



At a very large optical depth, which I denote by % = N d , I set the flux equal to zero, 
i.e., v™ dj = 0, and applying equation (38) on the last half of the last grid cell, denoted by 
(iV d - 1/2, JV d ), I obtain 

JV d -l_ N d / \ 
es cs m \ r>m,s m ^m,s I 

*(Ps)n a -i,2jN [ UN ^ " \ RN ^ k " • (43) 

Note that this boundary condition is different than equation (16). Finally, at the first and 
last energy grid-points, boundary conditions (17) become simply 

«™l+(«-l)JV M = M^v e+ (/_i)at m = ^i + (i_i)JV„ = V i?N e +(l-l)N„ = ■ (44) 

The system of algebraic equations (39)-(44) can be written in the general matrix form 

A?ut + B?ut +1 = C?v? +1/2 - G?> a ,i = l,...,N d -l 
TCi/2 + W + i/2 = F?uf-GT S ,i = l,...,N d , (45) 

where A™, B?, D™, and £™ are N.N,, x N e N^, diagonal matrices, C™ and are full 
N e N^ x N e N^ matrices, and G™' a and G™' s are N e N^ x 1 vectors. Equations (45) can then 
be solved recursively from i — N d to i — 1, using the relations 

< = + (46) 

^1+1/2 — ^i+l/2' U ™ 1 + ; (47) 

where 

VT = (F™ - DTU^-'ET (48) 

^ = ( F ™ - TO/aJ'^TOi/j + Gr) (49) 

0*1/2 = (c^-^wr 1 ^ (so) 

^+1/2 = (cr-^T)" 1 WCv J + GT?/ 2 )- (51) 

Note here a misprint in the elimination scheme of Milkey et al. (1975) as well as a difference 
with the above scheme that arises because of the presence of the antisymmetric term G™'"/ 2 - 

It is important to point out here that the method presented in this section is very efficient 
since it requires the solution of only four equations of simplicity equal to the minimum 
required for calculating the interaction of an external radiation field with an accretion disk. 
Therefore, the computational cost of this method is exactly equal to four times the cost 
required for solving the simplest problem of normal illumination of a plane-parallel slab. As 
an example, the solution of a problem with 20 grid points in optical depth, 10 grid points 
in angle, and 60 grid points in energy (i.e., similar to the one shown in Fig. 1) requires only 
2 CPU minutes on a 500 MHz Alpha processor. 
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Fig. 1. — The energy dependence of the vertical flux emerging from an illuminated accretion disk, decom- 
posed into the various expansion orders of the specific intensity. The angle of illumination is set to 9- ln = 80° 
with respect to the normal and the photon index of the illuminating spectrum to Y = 2.5. 



4. APPLICATION TO ILLUMINATED, COLD ACCRETION DISKS 

In this section, I use the method described in §3 for calculating the albedo of an illumi- 
nated, cold accretion disk. I consider concentric annuli of the disk, which I approximate by 
plane-parallel slabs. I take into account Compton scattering as well as bound-free absorption 
from a cold gas (Morrison & McCammon 1983) and neglect all other radiation processes. 
When the radiation field is weak, the ratio of the bound-free absorption coefficient to the 
scattering cross section is independent of the electron density and temperature. As a result, 
I can solve the radiative transfer problem in terms of the electron-scattering optical depth 
r cs without the need to consider the vertical disk structure. When the radiation field is 
strong, the illumination of the disk will affect both its ionization balance and the absorption 
coefficients, and hence the overall solution will depend explicitly on the vertical structure of 
the disk (see, e.g., Nayakshin & Kallman 2001). 

In general, the radiative transfer equation (7) is linear in the specific intensity. I can, 
therefore, write 

I(t, /i, 0, e) = I v (r, fi, 0, e) + I^r, fi, 0, e) , (52) 
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Fig. 2. — Relative contribution to the vertical flux of the various expansion orders of the 
specific intensity, for different angles of illumination. The vertical flux is calculated for a 
photon energy of 20 keV and a power-law spectrum of illumination with a photon index 
r = 2.5. 

where I v and I\ describe the radiation field due to viscous heating and illumination respec- 
tively, and solve for the two radiation fields independently. Note here that decomposition (52) 
is only formally valid, since the absorption and emission coefficients, as well as the scattering 
kernel may depend on the electron density and temperature that are determined by the total 
radiation field. The model parameters for the calculations include the angle of illumination, 
9\ with respect to the normal and the spectrum of the illuminating radiation, which I set to 
a power-law of photon index T, i.e., F(e) ~ e~ r+1 . 

Figure 1 shows the contribution to the vertical flux (Eq. [21]) of the various orders of the 
decomposition (8), for an illumination angle of 9[ = 80° and a power-law spectrum of photon 
index T = 2.5. Even though the term of zeroth-order provides the dominant contribution 
to the vertical flux, the contribution of higher-order terms is not negligible. This is shown 
in Figure 2, where the contribution of the high-order terms to the vertical flux is plotted 
for different angles of illumination. The total correction, caused by the obliqueness of the 
illumination, ranges between ~ 5 — 20%. Note here that the relative contributions of the 
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different orders of the specific intensity plotted in Figure 2 correspond to the case of zero 
electron temperature and increase for finite electron temperatures because of the asymmetry 
of the Compton terms in the scattering kernel. 

In the context of accretion onto galactic compact objects, the albedo rj of an accretion 
disk is usually defined in terms of the fraction of the illuminating flux that does not heat 
the disk gas (see, e.g., de Jonk et al. 1996). In this section, I first give the general expression 
for the albedo of an accretion disk that is illuminated obliquely and then evaluate it for the 
case of geometrically thin, optically thick, cold accretion disks. 

Starting from equation (52) and integrating the transfer equation for /; over photon 
energy and directional angle I obtain 



V • #i = - X Ji ~ n e a T I 4— — ) Ji , 



m. 



(53) 



where Jj and Hi are the energy-integrated, zeroth and first angular moments of /; and 

= /eJj(e)de 
K) ~ / Ji(e)de ' 

The rate Q h at which the illuminating radiation heats the gas is therefore 



(54) 
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(55) 



and the albedo can be written in terms of the volume integral of Q as 
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4 rr~, 
A*i Jo 

2 r 

'ThFj 



X + n e a T 
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X + n e a T 
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F 
T 

- 1 - e 



de / d/x/^(r es ,e,/x) , (56) 



where F is the energy-integrated flux of the illuminating radiation and r eSimax is the vertical 
height of the accretion disk. Expression (56) allows the study of the vertical profile of the 
heating of the accretion disk by the irradiating atmosphere as well as the calculation of its 
albedo, by solving only the zeroth-order transfer equation for I^(t ss , e, /i). 

Figure 3 shows the albedo of a cold accretion disk, calculated in the 2 — 20 keV energy 
range, for different power-law spectra and angles of illumination. As expected, for larger 
photon indices of the illuminating radiation, the fraction of low-energy photons, which are 
efficiently absorbed, is higher and hence the albedo of the disk is lower. At the same time, as 
the angle of illumination increases, the photons interact with the electrons in a shallower layer 
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Fig. 3. — The albedo of a cold accretion disk, calculated in the 2 — 20 keV energy range, for 
different power-law spectra and angles of illumination. 



of the accretion disk and therefore have a higher chance of escaping after one interaction, 
increasing the value of the albedo. The overall effect of our treatment of the obliqueness of 
irradiation is this systematic increase of the disk albedo with illumination angle, which can 
be up to a factor of ~ 2 larger that in the case of normal illumination. 



5. CONCLUSIONS 

In this paper, I studied the transfer of radiation in an accretion disk that is obliquely 
illuminated by an external source of radiation. I showed that the resulting transport problem 
can be decomposed exactly into four one-dimensional problems, which I solved using a variant 
of the Feautrier method. I then applied this method in calculating the albedos of cold 
accretion disks. 

The calculated values for the albedos are < 0.30, even for the softer spectra and larger 
illumination angles considered here. These values are small and cannot account, for exam- 
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pie, for the observed optical magnitudes of galactic low-mass X-ray binaries, which require 
albedos in excess of ~ 0.90 (see de Jong et al. 1996). Figure 3 shows that relying on a 
near-grazing illumination of the accretion disk is not enough to account for the observations. 
High ionization fractions at the surface layers of the disk, which would reduce the absorp- 
tion of photons, or even the existence of a highly ionized scattering wind above the accretion 
disk is probably required for the calculated albedos to reach the high values inferred from 
observations. 

In this study, I assumed for simplicity that all heavy elements in the accretion disk are 
neutral and, therefore their interaction with the illuminating photons is described by the 
bound-free opacities of Morrison & McCammon (1983). In reality, however, the heated skin 
of the accretion disk will be collisionally- and photo-ionized and its vertical ionization and 
thermal balance will need to be calculated self-consistently with the radiation field (e.g., 
Raymond 1993; Ko & Kallman 1994). Note, however, that the calculation of both the 
ionization balance and the radiative equilibrium depend only on the zeroth moment of the 
specific intensity (see eq. [18] and [55]) and, therefore, require the solution of only the zeroth- 
order transfer equation. As a result, the properties of the disk gas can be calculated exactly 
in a simple, one-dimensional configuration and the full angular dependence of the radiation 
field can then be calculated with prescribed gas properties. The results of a self-consistent 
caclulation of the radiation and gas properties will be reported elsewhere. 

I am grateful to G. Rybicki for bringing to my attention the possibility of decompos- 
ing a multi-dimensional transfer equation into a small number of one-dimensional equations 
and for carefully reading the manuscript. I also thank Feryal Ozel for many useful dis- 
cussions, especially on the implementation of the Feautrier method in problems with no 
forward-backward symmetry. This work was supported by a postdoctoral fellowship of the 
Smithsonian Institution and also, in part, by NASA. 
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